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Ab initio based polarizable force field generation and application to liquid 
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We extend the program potfit, which generates effective atomic interaction potentials from ab initio data, 
to electrostatic interactions and induced dipoles. The potential parametrization algorithm uses the Wolf 
direct, pairwise summation method with spherical truncation. The polarizability of oxygen atoms is modeled 
with the Tangney-Scandolo interatomic force field approach. Due to the Wolf summation, the computational 
effort in simulation scales linearly in the number of particles, despite the presence of electrostatic interactions. 
Thus, this model allows to perform large-scale molecular dynamics simulations of metal oxides with realistic 
potentials. Details of the implementation are given, and the generation of potentials for Si02 and MgO is 
demonstrated. The approach is validated by simulations of microstructural, thermodynamic and vibrational 
properties of liquid silica and magnesia. 
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I. INTRODUCTION 

As oxides are important in many technological appli- 
cations, they have been object of many experimental 
investigations 1-3 as well as of theoretical calculations. 4-8 
Despite continuous advances, ab initio methods cannot 
yet tackle the study of various material properties due to 
length and time scale restrictions. Classical effective po- 
tentials greatly increase the accessible system sizes and 
times; the quality of such simulations, however, depends 
crucially on the interaction models used. In this work, 
we describe a method to derive effective potentials for 
oxides in large-scale simulations. 

Both force field generation and simulation of oxide sys- 
tems are computationally much more demanding than 
those of metals or covalent materials due to long-range 
electrostatic interactions. There have been many at- 
tempts to create effective interaction potentials for oxides 
in the last decades. We combine two methods that have 
been applied successfully: First, in addition to charges we 
regard higher moments. It has been shown for silica 9,10 
that polarization effects of the oxygen atoms have to be 
taken into account and yield adequate results. The po- 
tential model of Tangney and Scandolo 9 (TS) determines 
the dipole moments by iteration to a self-consistent so- 
lution. This polarizable ion model yields significantly 
better results for many properties compared to other 
approaches. 10 Second, long-range interactions must be 
handled correctly and efficiently in Molecular Dynamics 
(MD) simulations. There are several methods with the 
one of Ewald 11 as the most famous. In a recent study for 
silica, 12 linear scaling in the number of particles could be 
achieved by using the Wolf summation method 13 without 
significant loss of accuracy. By optimizing the damp- 
ing and truncation of the long-ranged potential while 
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maintaining energy conservation, simulations can be per- 
formed at a comparatively small force-field cutoff radius 
of 8 A. 

The program potfit 14 ' 15 generates an effective atomic 
interaction force field with the force matching method, 16 
i.e., by adjusting the potential parameters to optimally 
reproduce a set of reference data computed in first- 
principles calculations. We have extended the program to 
electrostatic interactions and implemented the TS model 
with Wolf summation in the optimization algorithm. In 
this way, it is now possible to generate force fields for 
arbitrary metal oxides and use them in MD simulations, 
notably in the MD code IMD, 17,18 where the TS model 
and Wolf summation were already implemented. 12 The 
entire approach and the implementation in potfit are de- 
scribed in Sec. II. 

This potential generation method was applied to 
reparametrize the original TS force field for liquid sil- 
ica (SiCb). Because of its technological significance and 
its high abundance in nature, 19 silica has been thor- 
oughly investigated experimentally, 1,20 by using ab initio 
calculations 4,5,21 and in simulations with empirical in- 
teraction potentials. 6,9,22-24 Hence silica is an ideal test 
case for the new potential generation. Furthermore, a 
force field for liquid magnesia is generated. Magnesia 
(MgO) also is one of the most important metal oxides. 
Not only is it one of the simpler oxides (making it a fre- 
quently studied test system), it is also both ubiquitous 
and of technological importance. Hence there have been 
many experimental 2 ' 3,25,26 and theoretical 7, 8,27-29 stud- 
ies. By generating a force field for magnesia, we expand 
the combined use of the TS model and Wolf summation 
and demonstrate its power beyond silica. 

We determine basic microstructural, thermodynamic 
and vibrational properties for silica and magnesia in com- 
parison to experimental and first principles results and 
confirm the quality of this new force field approach by 
error estimations. The results are shown in Sec. Ill and 
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Sec. IV respectively. In Sec. V we sum up our results and 
give an outlook. 



II. FORCE FIELD 

A. Tangney-Scandolo potential model 

The TS potential contains two contributions: a short- 
range pair potential of Morsc-Strctch (MS) form, and a 
long-range part, which describes the electrostatic interac- 
tions between charges and induced dipolcs on the oxygen 
atoms. The MS interaction between an atom of type i 
and an atom of type j has the form 



u!f s = Di 



exp[7y(l 



1 ',; 
,.() 



2exp[f (1 



a) 

with Tij = \rij\, Tij = rj — n and the model parameters 
Dij, jij and r° , which have to be optimized. 

The dipole moments depend on the local electric field 
of the surrounding charges and dipoles. Hence a self- 
consistent iterative solution has to be found. In the TS 
approach, a dipole moment pf at position r\ in iteration 
step n consists of an induced part due to an electric field 
E(r{) and a short-range part pf R due to the short-range 
interactions between charges qi and qj . Following Rowley 
et al, this contribution is given by 



SR 



E 

a** 



fij ( r ij ) 



with 



fc=0 
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(3) 
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SR 
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where ai is the polarizability of atom i and E(ri) the 
electric field at position rj, which is determined by the 
dipole moments Pj in the previous iteration step. Con- 
sidering the interactions between charges U qq , between 
dipole moments C/ pp and between a charge and a dipole 
U pq , the total electrostatic contribution is given by 



jjEL = uqq + uPq + jjpp 



and the total interaction is 



IP 



U MS + C/ EL . 



(5) 



(6) 



B. Wolf summation 

The TS potential consists of pairwise interactions only. 
The electrostatic energies of a condensed system are 
then described by functions with r~™ dependence, n € 
{1,2,3}. For point charges (V -1 ) it is common to apply 
the Ewald method, where the total Coulomb energy of a 
set of N ions, 



E qq 



N N 
J- \ ^ \ ^ gtgj 

2 2^ 2^ r 

i=i j=i lJ 



(7) 



is decomposed in two terms E qq and E^ q by inserting a 
unity of the form 1 = crfc(ftr) + erf (fir) with the error 
function 



erf(ztr) := 



dt i 



(8) 



The splitting parameter k controls the distribution of 
energy contributions between the two terms. The short- 
ranged erfc term is summed up directly, while the smooth 
erf term is Fourier transformed and evaluated in recip- 
rocal space. This restricts the technique to periodic sys- 
tems. However, the main disadvantage is the scaling of 
the computational effort with the number of particles in 
the simulation box, which increases as 0(iV 3 / 2 ), 30 even 
for the optimal choice of k. 

Wolf et al. 13 designed a method with linear scaling 
properties O(N) for Coulomb interactions. By taking 
into account the physical properties of the system, the 
reciprocal-space term E^ q is disregarded. It can be writ- 
ten as 



fij (r^j) was introduced ad hoc to account for multipole 
effects of nearest neighbors and is a function of very short 
range, bij is the reciprocal of the length scale over which 
the short-range interaction comes into play, Cij deter- where 
mines amplitude and sign of this contribution to the in- 
duced moment. Together with the induced part, one ob- 
tains 
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(10) 



is the charge structure factor. V is the volume of the sim- 
ulation box. The charge structure factor is the Fourier 
transform of the charge-charge autocorrelation function. 
In liquid systems and largely also in solids, the charges 
screen each other like in a cold dense plasma. This means 
that for small magnitudes of the wave vectors k, the 
charge structure factor is also small. In a previous study, 
we showed 12 that the reciprocal-space term is indeed neg- 
ligible compared to the real-space part, provided that the 
splitting parameter k is chosen small enough. As, how- 
ever, a smaller k results in a larger real-space cutoff r c , 
the latter has to be chosen two or three times the size 
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of a typical short-range cutoff radius in metals. So a 
compromise has to be found for MD simulations: On the 
one hand k has to be small to allow disregarding the 
reciprocal-space term, on the other hand it should not 
be too small to avoid an oversized cutoff radius. 

In addition, a continuous and smooth cutoff of 
the remaining screened Coulomb potential -E qq (r.y) = 
Qiqj erfc(fiTy)r^ is adopted at r c by shifting the poten- 
tial so that it goes to zero smoothly in the first two deriva- 
tives at r = r c . We use the Wolf method for charges and 
its extension to dipolar interactions. For more informa- 
tion about the Wolf summation of dipolc contributions 
and a detailed analysis of the energy conservation in MD 
simulations, see Ref. 12. 



C. Wolf dipole error estimation 

To confirm the validity of the Wolf summation for the 
dipole contributions, it has to be analyzed whether the 
dipolar reciprocal-space term can also be neglected. This 
approximation is justified in a system with no significant 
long-range dipole-correlations like in ferroelectrics. The 
reciprocal-space term of the total energy of TV dipole mo- 
ments Pi can be written - analogously to Eq. (9) - as 

fc^O 1 1 

where e is the elementary charge and Q(k) the dipole 
structure factor 

1 N 

9i k ) '■= jw~2 p*®Pj cxp(ife • nj), (12) 

with the normalization factor 1/V Ne 2 . In each system 
of interest, the reciprocal-space term has to be negligible 
compared to the real-space contribution, which can be 
checked by evaluating the dipole structure factor. 



D. Implementation in potfit 

The programm potfit generates an effective atomic in- 
teraction force field solely from ab initio reference struc- 
tures. The potential parameters are optimized by match- 
ing the resulting forces, energies and stresses to according 
first-principles values with the force matching method. 
All reference structures used in this study were pre- 
pared with the plane wave code VASP. 31,32 For N m par- 
ticles, reference configuration m provides one energy e^, 
six components of the stress tensor l (I = 1,2, ..,6) 
and 3-/V m total force cartesian components n (n = 
1, 2, 3N m ) on N m atoms. The function 

Z = w e Z e + w s Z s + Zf (13) 



is minimized with respect to the charges, the parameters 
of the MS interaction Eq. (1) and the dipolar contribution 
Eqs. (2)-(4). Here 

M 

Z e =3^2 N m (e m -e° m ) 2 , 

m— 1 
- M 6 

^ = 2EE^(v-<if. ( 14 ) 

m=l 1 = 1 
M 3N m 

Z f = (fm,n - /m,n) 2 ' 

m — 1 n—1 

and e m , s m ./ and f m , n are the corresponding values cal- 
culated with the parametrized force field. w e and w s are 
certain weights to balance the different amount of avail- 
able data for each quantity. In the following we assume 
M reference structures that all consist of the same num- 
ber of particles (N m = N), but in principle, potfit can 
handle different numbers of particles for each reference 
structure. The root mean square (rms) errors, 

AF e = , AF S = J^i and AF f = \f^, 

e V 3MN ' V MN } V MN 

(15) 

are first indicators of the quality of the generated force 
field. Their magnitudes are independent of weighting fac- 
tors, number and sizes of reference structures. For the 
minimization of the potential parameters, a combination 
of a stochastic simulated annealing algorithm 33 and a 
conjugate-gradient-like deterministic algorithm 34 is used. 

We have extended the repertory of interactions in pot- 
fit by implementing a long-range electrostatic pair-force- 
routinc including the TS force field with Wolf summa- 
tion. Earlier we had also implemented it into the MD 
code IMD, 12 so that potential generation and simulation 
can now be used successively, potfit is parallelized us- 
ing the standard Message Passing Interface (MPI) 35 by 
distributing the reference configurations on the proces- 
sors. The whole method is not limited to diatomic metal 
oxides, but can treat any oxide. 

potfit now accepts the TS force field parameters as val- 
ues to be optimized. The MS potential Eq. (1) is defined 
by three parameters for each pair of interaction partners 
ij. Including the charges the polarizability a of the 
oxygen atom and the parameters bij and (Eqs. (2)- 
(4)), which only differ from zero in the case i ^ j, there 
are 14 parameters for a binary oxide. The requirement of 
charge neutrality reduces the number of free parameters 
by one. Thus, there are 13 parameters to optimize. 

III. APPLICATION TO SILICA 

A. Parametrization 

The new potential generation method has first been 
applied to liquid silica. Wc prepared five MD trajec- 
tories using the Wolf summed TS force field without 
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TABLE I. Force field parameters for silica, given in IMD units 
set eV, A and amu (hence charges are multiples of the elemen- 
tary charge). 



rcparametrization. 12 These constant-temperature runs 
were performed at five different temperatures between 
2000 and 4000 Kelvin. In addition, the volume was 
slightly lowered to prepare different pressure conditions. 
Then we took several snapshots following 10 ps of equi- 
libration. In this way, we prepared 47 liquid reference 
structures with on average 109 atoms, altogether 5123 
atoms. The reference structures show a pressure spec- 
trum from zero up to 15 GPa. They are used as in- 
put configurations for the first-principles plane wave code 
VASP, 31,32 where PAW pseudopotentials 36 and a gener- 
alized gradient approximation (GGA) of the exchange- 
correlation functional were used. With the local-density 
approximation (LDA) exchange-correlation functional, 
the well-known underestimation of the volume yielded 
clear deviations from experimental data. The approach 
with GGA and the same reference structures, however, 
did not overestimate the volume in MD simulations. Us- 
ing ultrasoft pseudopotentials to generate a reference 
daase had no noticeable influence on the results. 

The weights in potfit were chosen to w e = 0.1 and 
w s = 0.5, which is consistent with comparable opti- 
mization approaches. 6,9 The resulting parameters, how- 
ever, remained rather unaffected by modification of the 
weights. By contrast, the splitting parameter k had to be 
optimized as described in Sec. II B. Setting k = 0.02 A , 
a cutoff radius of only 8 A was found to be sufficient, 
which is small compared to other 6,12,37 long-range poten- 
tial approaches and results in an additional speedup in 
simulations. For comparison only, the procedure was re- 
peated with a cutoff radius of 10 A. The obtained results 
were quite similar. As the computational effort scales 
with r 3 ut , simulations with a cutoff radius of 8 A are 
about two times faster than the same settings with a 
cutoff radius of 10 A (see Fig. 5). It has to be mentioned 

that K = 0.02 A 1 describes a relatively weak damping. 
This identifies the good native screening ability of liquid 
silica. 

The final set of parameters is shown in Table I. The 
rms errors arc AF e = 0.1922, AF S = 0.0341 and AF f = 
1.6211. 
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FIG. 1. Equation of state of liquid silica at 3100 K for 
the new potential compared with experiment 20 and ab ini- 
tio calculations. 4 

B. Wolf dipole error estimation 

To analyze reciprocal-space contributions we simulated 
liquid silica with 4896 atoms at 3100 K and averaged 
several observables over the full simulation time of 200 ps. 
The total dipole moment isp = 3.47 x 10 - 28 Cm, which is 
small compared to a fully polarized system and thus can 
be taken fluctuation. Furthermore, we calculated 
the dipole structure factor and the resulting reciprocal- 
space contribution of the total energy, given in Eq. (11). 

For the selected damping of n — 0.02 A 1 we get 

1^ = 0.09^1/, (16) 

whereas the real-space contribution of the total energy 
per atom is 

—El ot = 7.20 eF. (17) 

Hence the reciprocal-space term is indeed insignificant 
and can be neglected. 

C. Results 

We validated the potential by determining thermody- 
namic, microstructural and vibrational properties of sil- 
ica. The goal is to cover different temperature and pres- 
sure scenarios in order to show the transferability of the 
new potential. The following simulations were all per- 
formed with the same initial configuration consisting of 
4896 atoms (1632 Si and 3264 O). 
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FIG. 2. Radial distribution functions for Si-Si, Si-O and O-O 
at 3000 K compared with ab initio calculations. 4 




80 100 120 140 160 180 



Si-O-Si angle (degrees) 

FIG. 3. Oxygen centered angle distribution in liquid silica at 
2370 K for the new potential compared with an ab initio MD 
study. 5 



In Fig. 1, we show the equation of state of liquid sil- 
ica at 3100 K calculated with the new potential. Pres- 
sures were obtained as averages along constant-volume 
MD runs of 10 ps following 10 ps of equilibration. For 
comparison, emperimental data 20 and first principles re- 
sults using VASP 4 are illustrated. The equation of state 
calculated with the new potential coincides with ab initio 
and experimental results. 
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FIG. 4. VDOS of amorphous silica at 300 K calculated with 
the new potential compared with an ab initio MD study. 21 a) 
partial VDOS for silicon atoms, b) partial VDOS for oxygen 
atoms, c) generalized VDOS. 



The radial distribution functions for Si-Si, Si-0 and 
0-0 at 3000 K with volume/Si0 2 V Si o 2 = 45.80 A 3 are 
evaluated for 100 snapshots taken out of an 100 ps MD 
run. The averaged curves are given in Fig. 2. Results 
calculated with the new potential are in accurate agree- 
ment with ab initio data. 4 The Si-Si and 0-0 curves are 
slightly shifted to larger distances. 

The Si-O-Si angle distribution was determined from 
300 MD runs at 2370 K and zero pressure. The averaged 
result is compared with an ab initio MD study 5 and is 
depicted in Fig. 3. As discussed in Ref. 5, ab initio MD 
tends to shift the angle distribution of liquid silica to 
slightly smaller angles with respect to the distribution 
generated by empirical MD simulations. 

The new force field was also applied to simulations 
of amorphous silica, although its parameters were only 
optimized with liquid reference structures between 2000 
and 4000 K. For this purpose, the liquid initial structure 
was cooled down to 300 K at an annealing rate of 0.01 
K/fs, which is recommended by Vollmayr et al. 38 and also 
used in Ref. 6. Using different annealing rates, however, 
had no significant impact on the results. The partial vi- 
brational density of states (VDOS) G si (E) (G (E)) for 
silicon (oxygen) was obtained by computing the Fourier 
transform of the timc-dcpcndcnt velocity-velocity auto- 
correlation function from a 100 ps MD trajectory with 
the software package nMoldyn. 39 The generalized VDOS 
G(E) is then calculated by 

G(E)= V ^G,(E) (18) 
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c(A) 


Si-O (A) 


Si-0-Si(°) 


Now potential 


5.15 


5.50 


1.65 


148.5 


Theory 


4.97 40 


5.39 40 


1.61 41 


145 41 



TABLE II. Lattice parameters, Si-0 bond length and Si-O-Si 
angle of a-quartz compared with theoretical studies. 

with the scattering cross section and atomic mass 
of atom fi. Fig. 4 shows the partial and general- 
ized VDOS compared with an ab initio study from Ref. 
21. The curves were adjusted in order to display the 
three main bands of each curve at similar frequencies. To 
achieve this, a global, constant relative frequency scaling 
for the MD curves, 

w = (1+7)0/ (19) 

was introduced, where 0/ is the original eigenfrequency, 
u! the scaled frequency, and 7 the constant relative fre- 
quency scaling. For silica we find the optimal 7 = 0.3. 
The new potential reproduces the key features of the 
partial VDOS. There are, however, several smaller devia- 
tions from ab initio results: The partial VDOS for silicon 
overestimates the low-energy band as the main-peak is 
shifted by around 12 meV to lower energies. The partial 
VDOS for ogyxen does not reproduce the band between 
35-55 meV and the two smaller peaks at around 77 meV 
and 108 meV. These characteristics are also reflected in 
the generalized VDOS. The relatively large scaling of the 
MD curves might also be a weakness of the new potential. 
In summary, the new potential is able to qualitatively re- 
flect the lattice dynamics, although it was not optimized 
for amorphous states at 300 K, but there are limits to its 
accuracy. 

As a final test, we used the force field for MD studies 
of a-quartz, one of the most important low-pressure crys- 
tal structures of Si02. This is a very hard test for the 
transferability of the potential, considering liquid struc- 
tures between 2000 and 4000 K were used as reference 
data. In Table II, we compare the lattice parameters, 
Si-0 bond length, and the Si-O-Si angle to analytical 
values. 40,41 The force field overestimates the equilibrium 
volume; the lattice parameters are too large by on av- 
erage 2.8%. In an MD simulation at 300 K over 10 ns, 
our force field stabilized the alpha-quartz structure. This 
shows that the new potential does yield reasonable results 
even under conditions it was not optimized for. 

We simulated systems with up to 2.5 million atoms 
in order to probe the linear scaling properties of the 
model in IMD. To judge the performance of the Wolf 
summation, the CPU cost of simulating a system of 
point charges was compared to the P3M 42 (Particle- 
particle/particle- mesh) method, as implemented in the 
MD code ESPResSo. 43 In P3M, the charges associated 
with the smooth erf term (see Sec. II B) are transferred 
to a lattice, where then the Poisson equation is solved. 
The CPU time per time step and atom is shown in Fig. 5. 
It can be seen clearly that the computational effort with 
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FIG. 5. Scaling of computational effort with system size for 
various long-range interaction techniques. For pure point 
charges, P3M is fastest for systems with less than 80 000 
atoms. In larger systems, Wolf summation performs better. 
The TS model also scales linear with system size using Wolf 
summation (it is not implemented in Ewald summation or the 
P3M method); it is slower by a factor of about 2.6 compared 
to point charges. 

Wolf summation scales perfectly linear in our implemen- 
tation. For smaller systems (fewer than about 80 000 
particles), ESPResSo is faster, but as the systems become 
larger, the 0(N log N) scaling of P3M loses to the linear 
scaling of Wolf. As a reference, we also show the CPU 
cost of the Ewald method in IMD. All calculations were 
performed on a single 2.83 GHz Intel Nehalcm CPU. The 
computational cost per atom of the TS model is indepen- 
dent of system size; compared to coulombic charges, the 
TS model in Wolf summation is slower by a factor of 2.6. 
This shows that in general, the number of steps in the 
self-consistency loop is independent of system size. 

IV. APPLICATION TO MAGNESIA 

A. Parametrization and error estimation 

To generate the reference structure database for liquid 
magnesia, we used snapshots from an MD trajectory with 
an ad-hoc potential. We prepared 80 liquid reference 
structures with on average 242 atoms, in total 19 360 
atoms. The reference forces, stresses and energies were 
computed with VASP using PAW pseudopotentials and 
in the GGA approximation. From these, we obtained an 
intermediate TS potential, which then was again used 
to generate a new MD trajectory and then a new set of 
reference structures. This procedure was iterated until 
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TABLE III. Force field parameters for magnesia, given in IMD 
units set eV, A and amu. 



there was no further significant change in the potential 
parameters. The pressure spectrum of the final reference 
structure database is between and 15 GPa, while the 
temperature varies between 2000 and 5000 K, to account 
for the higher melting point of magnesia. 

The weights in potfit again were set to w e = 0.1 and 
w s = 0.5. A cutoff radius of 8 A was sufficient when 
choosing k = 0.1 A . The final set of parameters is 
shown in Table III. The rms errors are AF e = 0.1170, 
AF S = 0.0228 and AF } = 0.5295. 

The following values for the Wolf dipole error estima- 
tion again were obtained by averaging over a simulation 
time of 200 ps. In liquid magnesia with 5832 atoms at 
5000 K, the total dipole moment is p = 4.76 x 10~ 33 
Cm, which is even smaller than that of silica. Also the 
reciprocal-space term is again really insignificant. For 
k = 0.1 A we get 



1 

N 



£PP = 0.15 fieV, 



(20) 



whereas the real-space contribution of the total energy 
per atom is 



1 



6.30 eV. 



(21) 



B. Results 

The new potential is validated by determining thermo- 
dynamic and microstructural properties of liquid mag- 
nesia. We also investigated the transferability of the 
new potential beyond the optimized temperature range 
by modelling the most important crystal structure per- 
iclase (NaCl-type). The following simulations were all 
performed with the same initial configuration consisting 
of 5832 atoms (2916 Mg and 2916 O). 

The equation of state, obtained as in silica, is shown 
in Fig. 6. Due to the high melting point of magnesia we 
chose a temperature of 5000 K. The result is compared 
with a first- principles study using VASP. 7 Although the 
new potential is optimized with a reference database hav- 
ing a pressure spectrum from zero to 15 GPa, it is able to 
reproduce the ab initio results quite accurate up to 160 
GPa. 
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FIG. 6. Equation of state of liquid magnesia at 5000 K for 
the new potential compared with ab initio calculations. 7 
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FIG. 7. Radial distribution functions for Mg-Mg (which is 
very similar to O-O case) and Mg-0 at 3000 K, compared 
with ab initio calculations. 7 



In Fig. 7 the radial distribution function for Mg-Mg 
(which is very similar to the 0-0 function) and Mg-0 is 

depicted at 3000 K with volume/MgO V Mg o = 27.76 A 3 . 
Our results coincide precisely with ab initio data. 7 
The Mg-O-Mg angle distribution was determined from 

S 3 

several MD runs at T^MgO = 33.99 A and three differ- 
ent temperatures 3000 K, 4000 K and 5000 K. The curves 
look similar, they show a maximum at 100°. So the inter- 
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FIG. 8. Oxygen centered angle distribution in liquid magnesia 
for the new potential for 3000 K, 4000 K and 5000 K. 
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FIG. 9. Generalized VDOS of periclase at 300 K calculated 
with the new potential compared with ab initio calculations 44 
and an experimental study. 45 



atomic angles in the magnesia melt are about 10° greater 
than in periclase. Both distribution studies follow from 
an averaging in the same way as for silica. 

Finally we applied the new force field to simulations of 
periclase at 300 K, although its parameters were only op- 
timized with liquid reference structures between 2000 and 
5000 K. The lattice constant is in good agreement with 
recent ab-inito and experimental studies, as presented in 



New Potential 4.214 



Exp.^ 


4.212 


Exp. a 


4.211 


ab initio (GGA) 8 


4.234 


ab initio (LDAf 


4.240 



TABLE IV. Lattice constant a [A] of periclase determined 
with the new potential compared to ab initio and experimen- 
tal studies. 



Table IV. 

In the same way as for silica, we obtained the general- 
ized VDOS of periclase. Fig. 9 shows a comparison with 
ab initio calculations 44 and an experimental 45 study. A 
frequency scaling of only 7 = 0.03 is applied. The new 
potential is able to qualitatively reproduce the key fea- 
tures, but there a two weak points: Firstly, the main 
peak at around 53 meV is shifted to lower frequencies by 
around 10 meV. Secondly, the peak at 36 meV - originat- 
ing from the partial VDOS for oxygen - is only weakly 
reproduced. In summary, the new force can give only 
qualitative results in this temperature range, for which 
it was not optimized. 

Not surprisingly, the scaling properties of simulations 
of magnesia are the same as those of silica: The com- 
putational effort scales linear in the number of particles. 
Also, the TS model is slower than constant charges by a 
factor independent of system size. 



V. CONCLUSION 

In this work, we presented the extension of the program 
potfit to electrostatic interactions. It can be used to gen- 
erate a force field for any condensed oxygen system by 
force matching to ab initio reference data. The system of 
interest only has to have an electrostatic screening dispo- 
sition to apply Wolf summation. We have demonstrated 
that the approach using the TS potential model yields 
accurate results for liquid silica and magnesia. Prob- 
ing solid phases, which were not included in the refer- 
ence database during the force field optimization, shows 
that qualitative studies are possible; quantitative analy- 
sis, however, should be done carefully. 

A new force field for crystalline a-alumina is currently 
in preparation 46 using the methods introduced in the 
present publication. In this process many crystalline 
structures are included in the reference database for op- 
timizing the potential parameters. This yields a highly 
accurate description of crystalline properties. 

The TS potential is a pure pair term potential. Hence 
simulations are less expensive compard to other ap- 
proaches with many-body potentials like a three-body 
interaction approach. Nevertheless the many-body char- 
acter of the electrostatics is indued in the self-consistent 
dipole iteration approach. Due to the Wolf summation, 
we achieve linear scaling in the number of particles, which 



9 



makes it possible to investigate larger system sizes and 
timescales. 

Generating a force field for liquid silica is an ideal 
test case for the new approach, because there are sev- 
eral recent theoretical studies 6,9,12 using the TS poten- 
tial model. Beyond silica, we have presented a new force 
field for liquid magnesia, which yields generally accurate 
results. Our work shows that a crude starting poten- 
tial suffices to determine high-quality potentials, which 
makes the new force field approach applicable for other 
systems of interest. 

In comparison to former simulation approaches with 
long-range interactions, we present force fields with mini- 
mal cutoff radius. This results in faster simulations. The 
short cutoff shows that liquid silica and magnesia can be 
adequately described by the Wolf approximation. 

We hope that our new force field generation approach 
enhances the investigation of condensed oxide systems, 
which could not yet be adequately investigated with MD 
simulations, due to large system sizes or boundary con- 
dition restrictions. 
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